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Abstract 

We define fDistances, which generalize Euclidean distances, squared distances, and 
log distances. The least squares loss function to fit fDistances to dissimilarity data is 
fStress. We give formulas and R/C code to compute partial derivatives of orders one 
to four of fStress, relying heavily on the use of Faa di Bruno’s chain rule formula for 
higher derivatives. 
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code. 


1 Introduction 

The multidimensional scaling (MDS) loss function fStress (Groenen, De Leeuw, and Mathar 
(1995)) is defined as 


<r{x) ■= J2Y, - fix'Ajx)) 2 ( 1 ) 

l<i<j<n 


for some real-valued /. The Wij are positive weights, the Sij are dissimilarities. This does not 
necessarily imply that they are actual dissimilarity judgments or measurements, they can be 
arbitrary transformations of such judgments or measurements as well. Both and 5^ are 
fixed and known numbers, the A t] are fixed and known symmetric matrices. Note there is no 
constraint that either the dissimilarities S tJ or the fDistances fix'A^x) are non-negative, or 
that / is increasing. Thus minimizing fStress can also be used, for example, to fit fDistances 
to similarities. 

Our notation is somewhat different from the standard MDS notation, so let’s make some 
translations. We use x = vec(X), with X the MDS configuration of n points in p dimensions. 
If we define u(i, s) = (s — 1) * n + i, then = x v ^ s y For most of our formulas it is easier 

to work with vectors than it is to work with matrices, so we will use x instead of X. It is 
also helpful for our software development, since matrices and higher-dimensional arrays are 
stored as vectors anyway. 

The rip x np matrices A tJ are defined using unit vectors e, and e 3 of length n, all zero except 
for one element that is equal to one. If you like, the e* are the columns of the identity matrix. 
We first define n x n matrices A t] by 

Aij • (e% ty ) (e, Gj ) , 

and use p copies of A l3 to make the direct sum 

A^ := A^ © • • • © Aij ■ 

v y 

p times 

Thus the squared Euclidean distance between point i and j in the configuration X is 
tr X'AijX = x'AijX. 

1.1 fDistances 

An fDistance is a function of the form g{x) = f(x'Ax ) for some real-valued /. Thus minimzinf 
fStress is fitting fDistances to dissimilarities, using weighted least squares. fStress properly 
generalizes the usual stress function initially proposed by Kruskal (1964a) and Kruskal 
(1964b), which uses g{x) = \/x'Ax. The sStress loss function of Takane, Young, and De 
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Leeuw (1977) uses the identity g(x) = x'Ax, and the IStress function of Ramsay (1977) and 
Ramsay (1982) uses the logarithm g(x) = \og{x'Ax). 

For g a general power, i.e. g(x) = ( x'Ax) r , we get rStress, sometimes also known as qStress, 
studied in a host of (unpublished, my fault) papers by De Leeuw, Groencn, and Pietersz 
(2006), Groenen and De Leeuw (2010), De Leeuw (2014), De Leeuw, Groencn, and Mail' 
(2016c), De Leeuw, Groenen, and Mair (2016a), De Leeuw, Groenen, and Mair (2016d), De 
Leeuw, Groenen, and Mair (2016b). LIRL’s and/or DOI’s are in the references section. 

In Groenen, De Leeuw, and Mathar (1995) first and second partials of fSt.ress are given. 
We add third and fourth order partials. It is unclear if the higher order partials have any 
practical applications. In De Leeuw, Groenen, and Mair (2016d) there are some applications 
of the second derivatives of rStress. We briefly mention some possible applications of our 
higher order partials to majorization (a.k.a. MM) algorithms. 

2 Partials, Partials Everywhere 

2.1 Derivatives of fDistances 

We give expressions for the first four derivatives of an arbitrary, but sufficiently many times 
differentiable, fDistance. In fact, we first address a more general problem, which in its turn 
is a special case of the first four terms of the multivariate Faa di Bruno formula (see, for 
instance, Constantine and Savits (1996), Leipnik and Pearce (2007)). We then apply those 
results to fDistances. 

2.2 Faa di Bruno 

Suppose h : lZ n —> 71, g : 7 Z" — > 77, and / : 77 — » 77, such that h(x) = f(g(x)). We will give 
expressions for the partials of h of orders up to four. Our formulas separate the parts that 
depend on / from the parts that depend only on g. 

2.2.1 First 

The first partials can be written in the form 

T>ih(x) = Vf(g(x))h u (x), 


with 

hii{x) := V. L g(x). 


2.2.2 Second 

Next 


Vijh(x) = Vf(g(x))h 2 i(x) + V 2 f(g(x))h 22 {x), 
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h- 2 l(x) := T>ijg{x), 
h 22 (x) := T> i g(x)'D j g(x). 


2.2.3 Third 

Next 

V ijk h(x) = Vf(g(x))h 3 i(x) + V 2 f(g(x))h 32 (x) + V 3 f(g(x))h 33 (x), 

with 

h 31 (x) := V ijk g(x), 

h 32 (x) := V j g(x)'D ik g(x) + V l g(x)V jk g(x) + V k g(x)V i:j g(x), 
h 33 (x) := V i g(x)V j g(x)Vg k (x). 

2.2.4 Fourth 

And finally 

V ijk ih(x) = Vf(g(x))h A1 {x) + V 2 f (g(x))h A2 (x) + V 3 f(g(x))h 43 {x) + V 4 f(g(x))hu(x), 
with 


K l(ar) := V ijkl g(x), 

h 42 (x) := T> j ig(x)'D ik g(x) + Vag(x)V jk g(x) 

+ T> j g(x)V ik ig(x) + T> i g(x)'D jk ig(x) + V k g(x)V ij ig(x) + V l g(x)'D ijk g(x), 
h 43 {x) := 'D u g(x)'D j g(x)Vg k (x) + V i g(x)V jl g(x)Vg k (x) + V i g(x)V j g(x)V kl g(x) 
+ V j g(x)V ik g{x)Vig(x) + V i g(x)'D jk g(x)T>ig(x) + V k g(x)T> ij g(x)T>ig(x), 
h M (x) = V i g{x)V j g(x)V k g{x)Vgi{x). 


2.2.5 Quadratic Forms 

If g(x) = x'Ax for some symmetric A, as it is in MDS, then T>ig(x) = 2 {Ax}i and 'D lJ g(x) = 
2 dij. Also V i j k g{x) = 0 and V ljk ig{x) = 0. Note that we write {Ax}i for element i of the 
vector Ax. Also note that the results in this section apply for any symmetric matrix A, even 
if it is not positive-semidefinite. 


For the first partials 


for the second partials 


Vjh(x) = 2 Vf(x'Ax){Ax}ii 


T>ijh(x ) = 2 T>f{x' Ax)a,ij + 4 V 2 f(x'Ax)){Ax}i{Ax}j, 
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and for the third partials 

T>ij k h(x) = 4 V 2 f(x'Ax){{Ax} i a jk + {Ax}ja ik + {Ax} k a i:j } + 8V 3 f(x'Ax){Ax} i {Ax} j {Ax} k . 
For the fourth partials we again write 

T> ijk ih(x) = T>f(g(x))h 41 (x) + V 2 f(g(x))h 42 (x) + V 3 f (g(x))h 43 (x) +V 4 f(g(x))h M (x), 
and now we have 


h 4 i(x) := 0, 

h'42 ) • 41 (lj](li k 4“ Q'il&jk }) 

hi 3 (x) := 8{aii{Ax}j{Ax} k + aji{Ax}i{Ax} k + a k i{Ax}i{Ax}j+ 
+ a ik {Ax}j{Ax}i + a jk {Ax}i{Ax}i + aij{Ax} k {Ax}i}, 
h M (x) := 16{Ax}i{Ax}j{Ax} k {Ax}i. 

2.3 Derivatives of Stress 

If we expand fStress we hnd, using notation used initially by De Leeuw (1977), 

cr(x ) — C — p(x) + v(x), 


with 

P( x ) '■= w ij 8i j f{x l Aijx), 

v(x) ■= \ Y.Y, W i j f 2 ( X 'A ij x), 

l<i<j»<n 

and with C a constant not depending on x. 

Clearly both p and g 2 are weighted sums of fDist-ances, where the fDistances in rj are the 
squares of the fDistances in p. Thus the partial derivatives are also weighted sums of the 
partial derivatives of the fDistances. Thus the previous formulas give us all we need. 

3 Implementation 

As in many of our previous recent reports, most of the computing is done in a C dialect 
conforming with the ,C() interface of R, but also fairly easible integrated into C programs 
that do not depend on the presence of R. 

At the moment we have implemented the following fDistances, which all are arbitrary powers 
of five base functions. 
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h{x) 

= {log (x'Ax)} p , 

(2) 

h{x) 

= {x'Ax} p , 

(3) 

h{x ) 

= {exp(a/kb;)} p , 

(4) 

h{x) 

j x 1 Ax 

(5) 

{ 1 + x' Ax J 

h{x) 

= {log(l + x 1 Ax)Y- 

(6) 


There are helper functions in C which evaluate the base functions and all four of their partials 
at a given point. There is also R glue for each of the helper functions, and the R function 
checkDQ uses the symbolic differentiation capabilities of R to check the helper functions. 

Partial derivatives of fDistances are evaluated using the multivariate Faa di Bruno results in 
the previous sections for the quadratic case with g(x) = x'Ax. The derivatives of the powers 
of the five base functions are computed from the derivatives of the base functions, using the C 
function f StressPower(), which implements a univariate version of Faa di Bruno. Since our 
formulas cover arbitrary powers of the five base functions they apply to both the fDistances 
in p and their squares in rj. 

The two basic C functions are f aa_di_bruno () and fStressFaaDiBrunoO. The first com¬ 
putes the first four partials of f {x'Ax) for general A, the second those of f{x'Ax) for the Re¬ 
used in MDS. In the second case we do not have to store A because there is a simple recipe to 
compute its elements whenever they are needed. Our C implementation is far from optimal, 
because we store full multidimensional arrays of partials without using their super-symmetry. 
We also fill the arrays with full loops over all indices, without using the sparseness of the Aij. 

The R function checkFO checks the faa_di_bruno() derivatives using numeric differen¬ 
tiation (Gilbert and Varadhan (2016)), and the R function faaR() provides the glue for 
faa_di_bruno. 

Finally there is the C function f StressPartials (), which makes the weighted sum of the 
partials of the fDistances. It has the R glue function fStressRO, and for the gradient and 
the Hessian it is checked against numerical derivatives by checks(). In our examples the 
first and second partials from fStressRO are the same as those from the numDeriv package. 
Since the higher derivatives of the base functions and their powers given the same results 
as those from deriv() and D(), we can be reasonable sure that the higher derivatives are 
also correct. Checking this, and possibly improving the efficiency of the code, are topics for 
further research (and for a future version of this paper). 

4 Discussion 

There is clearly some merit to the idea of fDistances, and to the general form of fStress. It 
covers many of the previous forms of least squares MDS. The first and second partials are 
not new. They can be used in gradient, Gauss-Newton, and Newton-Raphson minimization 
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algorithms and to draw pseudo-confidence ellipsoids (De Leeuw (2017)). The third and fourth 
order partials allow for a better local approximation of fStress, which at least theoretically 
can lead to faster algorithms and more precise confidence regions. 

In future research we will try to exploit convexity and concavity, and bounds on the derivatives, 
of fDistances and the implications for majorization algorithms to minimize fStress. 

5 Appendix: Code 

5.1 R Code 

5.1.1 fStress. R 

dyn.load( "fStress .so") 

fLogR <- function (x) { 
h <- 
• C( 

"fStressLog" , 
x = as.double (x), 
fO = as.double(O) , 
fl = as.double (0) , 
f2 = as.double (0) , 
f3 = as.double (0) , 
f4 = as.double (0) 

) 

return (list ( 
x = h$x, 
fO = h$f 0, 
fl = h$f 1, 
f 2 = h$f 2, 
f3 = h$f3, 
f4 = h$f 4 

)) 


fldeR <- function (x) { 
h <- 
• C( 

"f Stressldentity" , 
x = as.double (x), 
fO = as.double(O) , 
fl = as.double (0) , 
f2 = as.double (0) , 
f3 = as.double (0) , 
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f4 = as.double (0) 


return (list ( 
x = h$x, 
fO = h$f0, 
fl = h$f 1, 
f2 = h$f 2, 
f3 = h$f3, 
f4 = h$f 4 

)) 


fExpR <- function (x) { 
h <- 
• C( 

"fStressExponent" , 
x = as.double (x), 
fO = as.double(O) , 
fl = as.double (0) , 
f2 = as.double (0) , 
f3 = as.double (0) , 
f4 = as.double (0) 

) 

return (list ( 
x = h$x, 
fO = h$f 0, 
fl = h$f 1, 
f2 = h$f 2, 
f3 = h$f3, 
f4 = h$f 4 

)) 


fBndR <- function (x) { 
h <- 
• C( 

"f StressBounded" , 
x = as.double (x), 
fO = as.double(O) , 
fl = as.double (0) , 
f2 = as.double (0) , 
f3 = as.double (0) , 
f4 = as.double (0) 

) 


return (list ( 
x = h$x, 
fO = h$fO, 
fl = h$f 1, 
f2 = h$f 2, 
f3 = h$f3, 
f4 = h$f 4 

)) 


fLgoR <- function (x) { 
h <- 
• C( 

"f StressLogPlusOne" , 
x = as.double (x) , 
fO = as.double(O) , 
fl = as.double (0) , 
f2 = as.double (0) , 
f3 = as.double (0) , 
f4 = as.double (0) 

) 

return (list ( 
x = h$x, 
fO = h$f 0, 
fl = h$f 1, 
f2 = h$f 2, 
f3 = h$f3, 
f4 = h$f 4 

)) 


fPowR <- function (x, fNumber, pPower) { 
h <- 
• C( 

"fStressPower" , 

x = as.double (x), 

fNumber = as.integer (fNumber), 

ppower = as.double (pPower), 

fO = as.double(O) , 

fl = as.double (0) , 

f2 = as.double (0) , 

f3 = as.double (0) , 

f4 = as.double (0) 

) 
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return (list ( 
x = h$x, 
fO = h$fO, 
fl = h$f 1, 
f 2 = h$f 2 , 
f3 = h$f3, 
f4 = h$f4 

)) 

} 


faaR <- function(x, fNumber, pPower, a) { 
n <- length (x) 
h <- 
• C( 

"faa_di_bruno" , 

x = as.double (x), 

n = as.integer (n), 

fNumber = as.integer (fNumber - 1), 

pPower = as.double (pPower), 

a = as.double (a), 

ax = as.double (rep(0, n)), 

gx = as.double (0) , 

hO = as.double (0), 

hi = as.double(rep (0, n)), 

h2 = as.double (rep(0, n ~ 2)), 

h3 = as.double (rep (0, n ~ 3)), 

h4 = as.double (rep(0, n * 4)) 

) 

return (list ( 
gx = h$gx, 
hO = h$h0, 
hi = h$hl, 
h2 = h$h2, 
h3 = h$h3, 
h4 = h$h4 

)) 


fStressR <- function (x, w, delta, p, fNumber, pPower) { 
r <- length (x) 
m <- length (w) 
n <- round (r / p) 
hh <- 


• C( 
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"fStressPartials" , 

x = as.double (x), 

w = as.double (w) , 

delta = as.double (delta), 

n = as.integer (n), 

p = as.integer (p), 

fNumber = as.integer (fNumber - 1), 

pPower = as.double (pPower), 

stress = as.double (0), 

qdist = as.double (rep(0, m)), 

fdist = as.double (rep(0, m)), 

hi = as.double (rep(0, r)), 

h2 = as.double (rep(0, r ~ 2)), 

h3 = as.double (rep(0, r ~ 3)), 

h4 = as.double (rep(0, r ~ 4)) 

) 

return ( 
list ( 
x = x, 

stress = hh$stress, 
qdist = hh$qdist, 
fdist = hh$fdist, 
hi = hh$hl, 
h2 = hh$h2, 
h3 = hh$h3, 
h4 = hh$h4 

) 

) 

} 

5.1.2 check.R 

libraryC'numDeriv") 

fList <- list (function (x) 
log(x) , 
function (x) 
x, 

function (x) 
exp(x) , 
function (x) 
x / (1 + x) , 
function (x) 
log (1 + x)) 
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checkD <- function (expr, order, value) { 

DD <- function(expr , name, order = 1) { 
if (order < 1) 

stop(" ' order 1 must be >= 1") 
if (order == 1) 

D(expr, name) 
else 

DD(D(expr, name), name, order - 1) 

} 

dd <- DD(expr, "x", order) 

x <- value 

return (eval (dd)) 

} 

checkF <- function (x, fNumber, pPower, a) { 
f <- function (x, fNumber, pPower, a) { 
g <- sum (a * outer (x, x)) 
return ((fList[ [fNumber]] (g) ) ~ pPower) 

} 

hO <- f(x, fNumber, pPower, a) 
hi <- grad (f, 
x, 

fNumber = fNumber, 
pPower = pPower, 
a = a) 

h2 <- hessian (f, 
x, 

fNumber = fNumber, 
pPower = pPower, 
a = a) 

return (list ( 
x = x, 
hO = hO, 
hi = hi, 
h2 = h2 

)) 

} 

checks <- function (x, w, delta, p, fNumber, pPower) { 
f <- function (x, w, delta, p, fNumber, pPower) { 
r <- length (x) 
n <- round (r / p) 

d <- as.vector(dist(matrix(x, n, p))) ~ 2 
e <- fList[[fNumber]](d) ~ pPower 
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sum (w * (delta - e) ~ 2) / 2 


} 

hO <- 

f ( 

x, 

w = w, 

delta = delta, 

P = P> 

fNumber = fNumber, 
pPower = pPower 

) 

hi <- 
grad ( 

f, 
x, 

w = w, 

delta = delta, 

P = P> 

fNumber = fNumber, 
pPower = pPower 

) 

h2 <- 

hessian ( 

f, 
x, 

w = w, 

delta = delta, 

P = P> 

fNumber = fNumber, 
pPower = pPower 

) 

return (list ( 
x = x, 
hO = hO, 
hi = hi, 
h2 = h2 

)) 

} 

5.2 C Code 

5.2.1 fStress.h 


13 


#ifndef FSTRESS_H 
#define FSTRESS_H 

#include <math.h> 

#include <stdbool.h> 
#include <stdio.h> 

#include <stdlib.h> 


static 

inline 

int 

VINDEX (const 

int) 

static 

inline 

int 

MINDEX (const 

int , 

static 

inline 

int 

SINDEX (const 

int , 

static 

inline 

int 

TINDEX (const 

int , 

static 

inline 

int 

AINDEX (const 

int , 

static 

inline 

int 

CINDEX (const 

int , 


const 

int , 

const 

int) ; 


const 

int , 

const 

int) ; 


const 

int , 

const 

int) ; 


const 

int , 

const 

int, const 

int , 

const 

int , 

const 

int, const 

int , 


const int) ; 
const int, const ir 


static 

static 

static 

static 


inline double SQUARE (const double); 
inline double THIRD(const double); 
inline double FOURTH (const double) ; 
inline double FIFTH (const double) ; 


static 

static 

static 

static 

static 


inline double MAX(const double, const double); 
inline double MIN(const double, const double); 
inline int IMIN (const int, const int); 
inline int IMAX (const int, const int); 
inline int IMOD (const int, const int); 


static inline int ASEEK(const int, const int, const int, const int, const int, const ini 


static inline 

int VINDEX (const 

int 

i) 

{ return i - 

i; > 



static inline 

int MINDEX (const 

int 

i, 

const int j, 

const int 

n) 

{ 

return (i 

} 

- 1) + (j - 1) * 

n; 






static inline 

int AINDEX (const 

int 

i. 

const int j, 

const int 

k, 

const int n, 


const 

int 

m) 

{ 




return (i 

> 

- 1) + (j - 1) * 

n + 

(k 

- 1) * n * m: 

i 



static inline 

int CINDEX (const 

int 

i, 

const int j, 

const int 

k, 

const int 1, const int 


const 

int 

m, 

const int r) 

{ 



return (i 

- 1) + (j - 1) * 

n + 

(k 

- 1) * n * m 

+ (1 - 1) 

* 

n * m * r; 
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static inline int SINDEX (const int i, const int j, const int n) { 
return ((j - 1) * n) - (j * (j - 1) / 2) + (i - j) - 1; 

} 

static inline int TINDEX(const int i, const int j, const int n) { 

return ((j - 1) * n) - ((j - 1) * (j - 2) / 2) + (i - (j - 1)) - 1; 

} 

static inline double SQUARE (const double x) { return x * x; } 
static inline double THIRD(const double x) { return x * x * x; } 
static inline double FOURTH (const double x) { return x*x*x*x; ^ 
static inline double FIFTH(const double x) { return x*x*x*x*x; } 

static inline double MAX (const double x, const double y) { 
return (x > y) ? x : y; 

} 

static inline double MIN (const double x, const double y) { 
return (x < y) ? x : y; 

} 

static inline int IMAX (const int x, const int y) { return (x > y) ? x : y; } 

static inline int IMIN (const int x, const int y) { return (x < y) ? x : y; } 

static inline int IMOD (const int x, const int y) { 

return (((x % y) == 0) ? y : (x % y) ); 

} 


static inline int ASEEK(const int n, const int p, const int u, 
int value = 0; 
if (abs (i - j) < n) { 

for (int s = 1; s <= p; s++) { 

int k = (s - 1) *n; 

if ((i == (u + k)) && (j == (u + k))) { 
value = 1; 

} 

if ((i == (u + k)) && (j == (v + k))) { 
value = -1; 

} 

if ((i == (v + k)) && (j == (u + k))) { 
value = -1; 

} 


if ((i == (v + k)) && (j == (v + k))) { 


const int v, 


const int i, 
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value = 1; 

} 

} 

> 

return (value); 

> 


#endif /* FSTRESS_H */ 

5.2.2 fStress.c 
#include "fStress.h" 

void fStressLog(const double *x, double *f0, double *fl, double *f2, double *f3, 
double *f4) { 
double xx = *x; 

*f0 = log(xx); 

*fl = 1 / xx; 

*f2 = -1 / SQUARE(xx); 

*f3 = 2 / THIRD(xx); 

*f 4 = -6 / FOURTH(xx); 
return; 

> 

void fStressIdentity (const double *x, double *f 0, double *fl, double *f2, 

double *f3, double *f4) { 

double xx = *x; 

*f0 = xx; 

*fl = 1.0; 

*f2 = 0.0; 

*f 3 = 0.0; 

*f 4 = 0.0; 
return; 

> 

void fStressExponent (const double *x, double *f0, double *fl, double *f2, 

double *f3, double *f4) { 

double xx = *x; 

*f0 = exp(xx); 

*fl = exp(xx); 

*f 2 = exp(xx); 

*f 3 = exp(xx); 

*f4 = exp(xx); 
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return; 

} 

void fStressBounded (const double *x, double *f0, double *fl, double *f2, 

double *f3, double *f4) { 
double xx = *x, xz = 1.0 + xx; 

*f0 = xx / xz; 

*fl = 1.0 / SQUARE(xz); 

*f2 = -2.0 / THIRD(xz); 

*f3 = 6.0 / FOURTH(xz); 

*f 4 = -24.0 / FIFTH(xz); 
return; 

} 

void fStressLogPlusOne(const double *x, double *f 0, double *fl, double *f2, 

double *f3, double *f4) { 
double xx = *x, xp = 1.0 + xx, xz = 1 / xp; 

*f0 = log(xp); 

*fl = xz; 

*f2 = -SQUARE(xz); 

*f3 = 2.0 * THIRD(xz); 

*f4 = -6.0 * FOURTH(xz); 

} 

void (*fStressTable [5] ) (const double *, double *, double *, double *, double *, 

double *) = {fStressLog, fStressIdentity, 

fStressExponent, fStressBounded, 
f StressLogPlusOne}-; 

void fStressPower(const double *x, const int *fNumber, const double *ppar, 

double *g0, double *gl, double *g2, double *g3, double *g4) { 
double fO, fl, f2, f3, f4, pp = *ppar; 

(void) fStressTable[*fNumber](x, &f0, &fl, &f2, &f3, &f4); 

*g0 = pow(f0, pp); 

*gl = pp * fl * pow(f0, pp - 1.0); 

*g2 = pp * (pp - 1.0) * pow(f0, pp - 2.0) * SQUARE(fl); 

*g2 += pp * pow(f0, pp - 1.0) * f2; 

*g3 = pp * (pp - 1.0) * (pp - 2.0) * pow(f0, pp - 3.0) * THIRD(fl); 

*g3 += 3.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * fl * f2; 

*g3 += pp * pow(f0, pp - 1.0) * f3; 

*g4 = pp * (pp - 1.0) * (pp - 2.0) * (pp - 3.0) * pow(f0, pp - 4.0) * 

FOURTH(fl); 

*g4 += 6.0 * pp * (pp - 1.0) * (pp - 2.0) * pow(f0, pp - 3.0) * SQUARE(fl) * 
f 2; 
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*g4 += 4.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * (fl * f3); 

*g4 += 3.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * SQUARE(f2); 

*g4 += pp * pow(f0, pp - 1.0) * f4; 

} 

void fStressFaaDiBruno(const double *x, const int *n, const int *p, 

const int *u, const int *v, const int *fNumber, 
const double *par, double *ax, double *gx, double *h0, 
double *hl, double *h2, double *h3, double *h4) { 
double fO, fl, f2, f3, f4; 

int nn = *n, uu = *u, vv = *v, pp = *p, np = nn * pp; 

*gx = 0.0; 

for (int i = 1; i <= np; i++) { 
ax [VINDEX(i)] =0.0; 

> 

for (int s = 1; s <= pp; s++) { 

double xuv = x[MINDEX(uu, s, nn)] - x[MINDEX(vv, s, nn)]; 
ax [MINDEX(uu, s, nn)] = xuv; 
ax [MINDEX(vv, s, nn)] = -xuv; 

> 

for (int i = 1; i <= np; i++) { 

*gx += ax[VINDEX(i)] * x[VINDEX(i)]; 

> 

(void)fStressPower(gx, fNumber, par, &f0, &fl, &f2, &f3, &f4); 

*h0 = fO; 

for (int i = 1; i <= np; i++) { 

hi [VINDEX(i)] = 2.0 * fl * ax[VINDEX(i)]; 

> 

for (int i = 1; i <= np; i++) { 

for (int j = 1; j <= np; j++) { 

h2[MINDEX(i, j, np)] = 2.0 * fl * ASEEK(nn, pp, uu, vv, i, j) + 

4.0 * f2 * ax[VINDEX(i)] * ax[VINDEX(j)]; 

} 

> 

for (int i = 1; i <= np; i++) { 

for (int j = 1; j <= np; j++) { 

for (int k = 1; k <= np; k++) { 
h3[AINDEX(i, j, k, np, np)] = 

4.0 * f2 * 


((ax [VINDEX(i)] 

* ASEEK(nn, 

PP> 

uu, 

vv, 

j. 

k)) + 

(ax[VINDEX(j)] 

* ASEEK(nn, 

PP> 

uu, 

vv, 

i, 

k)) + 

(ax [VINDEX(k)] 

* ASEEK(nn. 

PP> 

uu, 

vv, 

i. 

j))); 


h3[AINDEX(i, j, k, nn, nn)] += 

8.0 * f3 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax[VINDEX(k)]; 
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> 


} 

} 

> 

for (int i = 1; i <= np; i++) { 

for (int j = 1; j <= np; j++) { 

for (int k = 1; k <= np; k++) { 

for (int 1 = 1; 1 <= np; 1++) { 

h4[CINDEX(i, j, k, 1, np, np, np)] = 
4.0 * f2 * 


(ASEEK(nn, pp. 

uu, 

vv, 

j > 

1) * 

ASEEK(nn, 

PP> 

uu. 

vv. 

. i, 

ASEEK(nn, pp. 

uu, 

vv. 

i. 

1) * 

ASEEK(nn, 

PP> 

uu, 

vv. 

. j> 


— , j j -- y — ? j j ^ -* 

16.0 * f4 * ax[VINDEX(i)] * ax[VINDEX(j)] 
ax[VINDEX(k)] * ax[VINDEX(l)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, mi, vv, i, 1) * 
ax[VINDEX(k)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, uu, vv, j, 1) * 
ax[VINDEX(k)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, uu, vv, k, 1) * 
ax[VINDEX(j)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, uu, vv, i, k) * 
ax[VINDEX(l)] ; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, uu, vv, j, k) * 
ax[VINDEX(l)] ; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * ASEEK(nn, pp, uu, vv, i, j) * ax[VINDEX(k)] * 
ax[VINDEX(l)]; 


ax [VINDEX(j)] * 

ax [VINDEX(i)] * 

ax [VINDEX(j)] * 
ax [VINDEX(j)] * 

ax [VINDEX(i)] * 


return; 


void faa_di_bruno( const double *x, const int *n, const int *fNumber, 

const double *par, const double *a, double *ax, double *gx, 
double *h0, double *hl, double *h2, double *h3, double *h4) { 
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double fO, fl, f2, f3, f4; 
int nn = *n; 

*gx = 0.0; 

for (int i = 1; i <= nn; i++) { 
ax [VINDEX(i)] =0.0; 
for (int j = 1; j <= nn; j++) { 

ax[VINDEX(i)] += a[MINDEX(i, j, nn)] * x[VINDEX(j)]; 

} 

*gx += ax[VINDEX(i)] * x [VINDEX(i)]; 

> 

(void)fStressPower(gx, fNumber, par, &f0, &fl, &f2, &f3, &f4); 

*h0 = fO; 

for (int i = 1; i <= nn; i++) { 

hi[VINDEX(i)] = 2.0 * fl * ax[VINDEX(i)]; 

> 

for (int i = 1; i <= nn; i++) { 

for (int j = 1; j <= nn; j++) { 

h2[MINDEX(i, j, nn)] = 2.0 * fl * a[MINDEX(i, j, nn)] + 

4.0 * f2 * ax[VINDEX(i)] * ax[VINDEX(j)]; 

} 

> 

for (int i = 1; i <= nn; i++) { 

for (int j = 1; j <= nn; j++) { 

for (int k = 1; k <= nn; k++) { 
h3[AINDEX(i, j, k, nn, nn)] = 

4.0 * f2 * 

((ax[VINDEX(i)] * a[MINDEX(j, k, nn)]) + 

(ax[VINDEX(j)] * a[MINDEX(i, k, nn)]) + 

(ax[VINDEX(k)] * a[MINDEX(i, j, nn)])); 
h3[AINDEX(i, j, k, nn, nn)] += 

8.0 * f3 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax [VINDEX(k)]; 

} 

} 

} 

for (int i = 1; i <= nn; i++) { 

for (int j = 1; j <= nn; j++) { 

for (int k = 1; k <= nn; k++) { 

for (int 1 = 1; 1 <= nn; 1++) { 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] = 

4.0 * f2 * 

(a[MINDEX(j, 1, nn)] * a[MINDEX(i, k, nn)] + 
a[MINDEX(i, 1, nn)] * a[MINDEX(j, k, nn)]); 
h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

16.0 * f4 * ax [VINDEX(i)] * ax[VINDEX(j)] * 
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} 

} 

> 

return; 


ax[VINDEX(k)] * ax[VINDEX(l)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] 

8.0 * f3 * a[MINDEX(i, 1, nn)] * ax[VINDEX(j)] * 
ax[VINDEX(k)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * a[MINDEX(j, 1, nn)] * ax[VINDEX(i)] * 
ax[VINDEX(k)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * a[MINDEX(k, 1, nn)] * ax[VINDEX(j)] * 
ax[VINDEX(j)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * a[MINDEX(i, k, nn)] * ax[VINDEX(j)] * 
ax[VINDEX(l)]; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * a[MINDEX(j, k, nn)] * ax[VINDEX(i)] * 
ax[VINDEX(l)] ; 

h4[CINDEX(i, j, k, 1, nn, nn, nn)] += 

8.0 * f3 * a[MINDEX(i, j, nn)] * ax[VINDEX(k)] * 
ax[VINDEX(l)]; 


+= 


void fStressPartials(const double *x, const double *w, const double *delta, 

const int *n, const int *p, const int *fNumber, 
const double *par, double *stress, double *d, double *fd, 
double *fl, double *f2, double *f3, double *f4) { 
double par2 = 2 * *par, nn = *n, pp = *p, np = nn * pp, gx, hx, hO, gO; 
double *ax = (double *)calloc((size_t)np, sizeof(double) ); 
double *hl = (double *)calloc((size_t)np, sizeof(double) ); 
double *h2 = (double *)calloc( (size_t) SQUARE(np), sizeof(double) ); 
double *h3 = (double *)calloc((size_t)THIRD(np), sizeof(double) ); 
double *h4 = (double *)calloc( (size_t) FOURTH(np), sizeof(double) ); 
double *gl = (double *)calloc((size_t)np, sizeof(double) ); 
double *g2 = (double *)calloc( (size_t) SQUARE(np), sizeof(double) ); 
double *g3 = (double *)calloc((size_t)THIRD(np), sizeof(double) ); 
double *g4 = (double *)calloc( (size_t) FOURTH(np), sizeof(double) ); 

^stress = 0.0; 

for (int j = 1; j <= nn - 1; j++) { 

for (int i = j +1; i <= nn; i++) { 
int k = SINDEX(i, j, nn); 
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> 


(void)fStressFaaDiBruno(x, n, p, &i, &j , fNumber, par, ax, &hx, fehO, 

hi, h2, h3, h4); 

d [k] = hx; 
fd[k] = hO; 

^stress += w[k] * SQUARE (delta [k] - fd[k]); 

(void)fStressFaaDiBruno(x, n, p, &i, &j, fNumber, &par2, ax, &gx, 

&g0, gl, g2, g3, g4); 
for (int r = 1; r <= np; r++) { 
int ind = VINDEX(r); 

fl[ind] += w[k] * (0.5 * gl [ind] - delta[k] * hi [ind]) ; 

for (int s = 1; s <= np; s++) { 

int ind = MINDEX(r, s, np); 

f 2 [ind] += w[k] * (0.5 * g2[ind] - delta [k] * h2 [ind]); 
for (int t = 1 ; t <= np; t++) { 

int ind = AINDEX(r, s, t, np, np); 

f3[ind] += w[k] * (0.5 * g3 [ind] - delta[k] * h3[ind]); 

for (int u = 1; u <= np; u++) { 

int ind = CINDEX(r, s, t, u, np, np, np); 
f4[ind] += 

w[k] * (0.5 * g4[ind] - delta[k] * h4[ind]); 

} 

> 

> 


} 

} 

> 

^stress /= 2.0; 
free(ax); 
free(hl); 
free(h2); 
free(h3); 
free(h4); 
free(gl); 
free(g2); 
free(g3); 
free(g4); 
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